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Abstract 

, Hamiltonian Boundary Value Methods are a new class of energy preserving one step methods for 

the solution of polynomial Hamiltonian dynamical systems. They can be thought of as a generalization 
of collocation methods in that they may be defined by imposing a suitable set of extended collocation 
conditions. In particular, in the way they are described in this note, they are related to Gauss 
collocation methods with the difference that they are able to precisely conserve the Hamiltonian 
pH function in the case where this is a polynomial of any high degree in the momenta and in the 

generalized coordinates. A description of these new formulas is followed by a few test problems 
showing how, in many relevant situations, the precise conservation of the Hamiltonian is crucial to 
simulate on a computer the correct behavior of the theoretical solutions. 



1 Introduction 

(N . 

, Hamiltonian Boundary Value Methods (HBVMs) form a subclass of Boundary Value Methods (BVMs) , 

whose main feature is that of precisely conserving the Hamiltonian function associated with a canonical 
Hamiltonian system 

f y = JVH(y), 



(N 
O 
O 



X 




1 y(t ) = yo e R 2m , 1 ' " ' 

(J is the identity matrix of dimension m), in the case where such function is of polynomial type. 

Two key ideas have permitted the realization of HBVMs: the definition of discrete line integral and 
what we called extended collocation conditions. The former, first introduced in [151 116], represents the 
discrete counterpart of the line integral defined over conservative vector fields, while the second is a 
relaxation of the classical collocation conditions which assures the conservation of the energy along the 
numerical solution {y n } generated by the method itself. 

Just as an initial clarification, we briefly show how this new approach to the problem reads when the 
classical Gauss collocation method is considered (see [TBI Remark 2.1] for more details). Given a stepsize 
h > and a set of s abscissae ci < • • • < c s disposed according to a Gauss-Legendre distribution on [0, 1], 
the Gauss method of order 2s is defined by means of the following polynomial collocation problem: 

cr(to) = yo, 

(2) 

a (t + Cih) = J\7H(cr(t + c t h)), i = 1, . . . , s. 
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As is well known, conditions @ uniquely define a polynomial a(t) of degree s which is used to advance 
the solution by posing y\ = a (to + h), while the internal stages satisfy Yi = a(to + Cih), i = 1, . . . s. The 
coefficients of the Butcher array and the weights are given by 

bj = / £j(c)dc, ciij = / £j(c)dc, with £j(c) = TT — . 

JO rit j '': ~ C r 

The s-degree polynomial a(t) may be thought of as a path in the phase space linking the state vectors 
yo to y\ and passing through the stages {Yi}. Due to the conservative nature of the vector field, we have 
that 

H( yi ) ~ H(y ) = f VH(y) ■ dy = h [ &(t + rh) T WH(a(t a + T h))dr. (3) 

J a JO 

Now, the above integral is exactly computed by the Gauss quadrature formula with abscissae {a} and 
weights {hi} if the degree of the integrand is not greater than 2s — 1 which means that the degree of 
H(y), say v, must not exceed 2 (linear or quadratic Hamiltonians only). Under this assumption, taking 
into account the collocation conditions ([2]), we obtain 

s s 

H( Vl ) - H(y Q ) = hY,h(cT(ti)) T VH( n (U)) = -h^ biV T H(a(ti)) JVH(^(U)) = 0, (4) 

i=l i=l 

where ij = to + c %h- Thus, by following a different route, we have obtained the classical result that the 
Gauss methods conserve quadratic Hamiltonian functions while fails to conserve polynomial Hamiltonian 
functions of higher degree0 

The above example is the starting point of our approach: the discrete line integral is the first sum in 
((4]), which turns out to vanish for quadratic Hamiltonians, due to the collocation conditions ([2]). 

The next section reports a descriptive introduction to HBVMs with much emphasis to the key ideas 
they rely on. We refer the reader to the papers [31 |H HH1 [21 G3 EH HH Q] for the details about the basic 
theory and implementation of HBVMs, and to the monograph 6 a as a reference for the theory of BVMs. 

In Section [3] we report a number of test problems of some relevance in the literature, for which the 
precise conservation of the energy turns out to be a crucial feature for the correct reproduction of the 
long time behavior of the solutions. This will be testified by comparing HBVMs to Gauss methods which, 
by the way, are symplectic integrators. 

2 Hamiltonian Boundary Value Methods 

In this section we introduce HBVMs by slightly elaborating the arguments in [3J [H [S] . As said above, the 
basic idea which HBVMs rely on is the so called discrete line integral, which is the discrete counterpart of 
the line integral associated with a conservative vector field. In more detail, starting from ([3j , we consider 
a polynomial, of degree at most s, such that 

<r(to)=yo, (?(to + h)=yi, (5) 

providing an approximation to the solution on the interval [to, to + h]. We consider the following expan- 
sions, 



a (t + ch) = y^-Pj(c)7j', a(t Q + ch) = y Q + h V\ 3 / P 3 -(r)dr, 

Jo 



(6) 

4=1 J'=l 
where the (vector) coefficients {ji} are to be determined. We also assume that the polynomials {Pi} 
constitute an orthonormal basis, on the interval [0, 1], for the vector space n s _i of polynomials of degree 



1 This argument may be generalized to other classes of collocation methods. 
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at most s — 1, i.e., 

»i 



(i 



P i( T ) P j( T ) = h3 = 1, • • • ) S, 



with Sij the Kronecker symbol. Such polynomials can be easily obtained by a suitable scaling of the 
shifted Legendre polynomials 5 . Substituting the first expansion in ([6]) into the line integral in (|3]), and 
requiring the resulting expression to vanish, then gives 



s „i 

$>J / P j(T)VH(a(t + rh))dr = 0, 



which is certainly satisfied by choosing 

7i= / P j {r)J^H{a{t +rh))dr, j = l,...,s. (7) 
Jo 

Multiplication of ([7]) by h J Q C Pj(x)dx and summation over j then gives, by virtue of the second expansion 
in ©, 

s pC pi 

a{t + ch) = y + ftV / P 3 {x)dx Pj{T)JVH{a{t + tK))At, c€[0,1], (8) 

j=1 JO JQ 

Let us now assume that H (y) is a polynomial of degree v. Consequently, the integral appearing at the 
right-hand side in ((SJ can be exactly discretized by a Gaussian formula over k Gauss-Legendre abscissae 
{ci}, which we shall consider hereafter, provided that 

k > y. (9) 
Let us denote by {u>i} the weights of the quadrature formula in the interval [0, 1], and set 

Hi = a(t + Cih), atj = Pj{x)dx, i = l,...,k, j = l,...,s. (10) 

Jo 

Consequently, ([H]) can be (exactly) discretized as: 

s k 

yi = Vo + hy j ajj }] uJjPj {ct)JVH(yt), i = l,...,k. (11) 
j=i i=i 

Definition 2.1 The set of equations ill}) , to be solved for the unknowns {yi}, defines an HBVM with k 
steps and degree s, in short HBVM(k,s). 

For such a method, the following properties hold true [1]: 

• it has order 2 s for all k > s; 

• it is symmetric and perfectly ^4-stable (i.e., its stability region coincides with the left-half complex 
plane, C~ [5]); 

• for fc = s, it reduces to the Gauss-Legendre method of order 2s; 

• it exactly preserves polynomial Hamiltonian functions of degree v, provided that ^ holds true. 

Remark 2.2 The actual implementation of HBVM(k, s) can be seen to result in the solution of a system 
of (block) size s, whatever is the value of k considered Consequently, if needed, large values of k 
can be easily considered. 
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The arguments in the previous remark, allow us to consider the limit formula of (|10|) - (|lip . in the case 
where H(y) is non-polynomial, as k — > oo. Clearly such a limit is given by formula (J8J), which, according 
to [1], is named HBVM(oo, s) or oo-HBVM of degree s. 

However, we emphasize that formula ((5J) becomes an operative method only after that a suitable 
discretization of the inner integral is considered and, replacing the integral by a quadrature formula with 
k nodes, leads back to a HBVM(fc, s) method. 

One can easily argue that, since in the non polynomial case the quadrature formula can approximate 
the corresponding integral with an arbitrary accuracy, under suitable regularity assumptions for H(y), a 
practical conservation of the energy may be obtained [HE]. The term "practical" means that, in many 
general situations, when k is high enough, the method makes no distinction between the function H(y) 
and its polynomial approximation, being the latter in a neighborhood of size e of the former, where e 
denotes the machine precision. 

We end this section by observing that, by differentiating both members of (0, one obtains 

s »1 

&(t + ch) = y2 p j( c ) / Pj(r)JVH{a(t + Th))dT, ce[0,l], 
which at the points {cj} provides, assuming H(y) to be a polynomial and k large enough: 

s -i 

&(t + c t h) = y2,Pj{ci) \ P;j(T)J\7H(<T(t +Th))dT, i = l,,..,k. 

Such formulae (the former being the limit of the latter as k — > oo) can be regarded as a kind of extended 
collocation conditions that generalize conditions @, according to [181 Section 2] (see also j3]). 



3 Numerical tests 

We present a few numerical test highlighting the good behavior of HBVMs in the long-time simulation of 
Hamiltonian systems. A direct comparison of HBVMs with Gauss methods is reported in order to better 
emphasize the stability properties of the former methods even when compared to a well known class of 
symplcctic formulae^] 

The use of a large stepsize of integration is a prerogative in long-time simulation of an evolutionary 
problem but, in general, one is forced to reduce h under a critical threshold in order to guarantee the 
qualitative behavior of the theoretical solution to be well reproduced by the numerical solution. From 
this point of view, we show that HBVMs allow the use of larger stepsizes than Gauss methods, which 
states that the conservation of the Hamiltonian function plays an important role in detecting the correct 
topological features of the solutions. 



3.1 Sitnikov's problem 

One of the main problems in Celestial Mechanics is to describe the motion of N point particles of positive 
mass {mi} moving under Newton's law of gravitation when we know their positions {qi} and momenta 
{pi} at a given time. Such a dynamical system, called the iV-body problem, is in the form (|T|), with 
Hamiltonian 

1 N II 112 N *-l 

i=l 1 i=l j=l 

2 As was seen in the previous section, the choice of Gauss methods has also been dictated by the fact that they represent 
the generating formulae of HBVMs when we use a Gauss distribution of the abscissae, namely the Gauss method of order 
2s coincides with HBVM(s,s). 
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with G the gravitational constant. While the two-body problem is completely solved in the sense that we 
can describe explicitly all its solutions (see, e.g., [12]), this is no more the case, for N > 3. Consequently, 
numerical simulation is of interest, in such a case. 

The Sitnikov problem is a particular configuration of the 3-body dynamics. In this problem two bodies 
of equal mass (primaries) revolve about their center of mass, here placed at the origin, in elliptic orbits in 
the (x, j/)-plane. A third, and much smaller body (planetoid), is placed on the z-axis with initial velocity 
parallel to this axis as well. 

The third body is small enough that the two body dynamics of the primaries is not destroyed. Then, 
the motion of the third body will be restricted to the z-axis and oscillating around the origin but not 
necessarily periodic. In fact, this problem has been shown to exhibit a chaotic behavior when the 
eccentricity of the orbits of the primaries exceeds a critical value that, for the data set we have used, is 
e ~ 0.725 (see Figure [TJ). 




Figure 1: The left picture displays the configuration of 3-bodies in the Sitnikov problem. To an eccentricity 
of the orbits of the primaries e = 0.75, there correspond bounded chaotic oscillations of the planetoid as 
is argued by looking at the space-time diagram in the right picture. 




Figure 2: Left picture: relative error \H(y n ) — H(yo)\/\H(yo)\ of the Hamiltonian function evaluated along 
the numerical solution of the HBVM(18,2) and the Gauss method. Right picture: relative error \M(y n ) — 
M(yo)\/\M(yo)\ of the angular momentum evaluated along the numerical solution of the HBVM(18,2) 
and the Gauss method. 

We have solved the Kepler problem with Hamiltonian function (Tl2|) by the Gauss method of order 4 
(HBVM(2,2)) and by HBVM(18,2) (order 4 and 18 steps), with the following set of parameters: 
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TV G mi m-i ms e d h tmax 



3 111 10~ 5 0.75 5 0.5 1500 

where e is the eccentricity, d is the distance of the apocentres of the primaries (points at which the two 
bodies are the furthest), h is the time-step and [0, imax] is the time integration interval. The eccentricity 
e and the distance d may be used to define the initial condition y = [q ,p ] (see [TH] for the details): 

<z = [-§, 0, 0, §, 0, 0, 0, o, io- 9 ], 

Po = [o, - ±Vw, o, o, iVio, o, o, o, ±]. 




_2 5I I I I I I I iqI I I I I I I 1 I 

-3 -2 -1 1 2 3 50 100 150 200 250 300 350 

Figure 3: The Sitnikov problem solved by the Gauss method of order 4, with stepsize h = 0.5, in 
the time interval [0,1500]. The trajectories of the primaries in the (x, y)-plane (left picture) exhibit a 
very irregular behavior which causes the planetoid to eventually leave the system, as illustrated by the 
space-time diagram in the right picture. 
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Figure 4: The Sitnikov problem solved by the HBVM(18,2) method (order 4), with stepsize h = 0.5, 
in the time interval [0,1500]. Left picture: the trajectories of the primaries are ellipse shape. The 
discretization introduces a fictitious uniform rotation of the (x, y)-plane which, however, does not alter 
the global symmetry of the system. Right picture: the space-time diagram of the planetoid on the z-axis 
displayed (for clearness) on the time interval [0,350] shows that, although a large value of the stepsize h 
has been used, the overall behavior of the dynamics is well reproduced (compare with the right picture 
of Figure [1]) . 

First of all, we consider the two pictures in Figure [5] reporting the relative errors in the Hamiltonian 
function and in the angular momentum evaluated along the numerical solutions computed by the two 
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Figure 5: Distance between the two primaries as a function of the time, related to the numerical solutions 
generated by the Gauss method (left picture) and HBVM(18,2) (right picture). The maxima correspond 
to the distance of apocentres. These are conserved by HBVM(18,2) while the Gauss method introduces 
patchy oscillations that destroy the overall symmetry of the system. 

methods. According to ©, we know that the HBVM(18,2) precisely conserves Hamiltonian polynomial 
functions of degree at most 18. This accuracy is high enough to guarantee that the nonlinear Hamiltonian 
function (|12[) is as well conserved up to the machine precision (see the left picture): from a geometrical 
point of view, this means that a local approximation of the level curves of (|12[) by a polynomial of degree 
18 leads to a negligible error. The Gauss method exhibits a certain error in the Hamiltonian function 
while, being this formula symplectic, it precisely conserves the angular momentum, as is confirmed by 
looking at the right picture of Figure [2l From the same picture, one sees that the error in the numerical 
angular momentum associated with the HBVM(18,2) undergoes some bounded periodic-like oscillations. 

Figures [3] and 0] show the numerical solution computed by the Gauss method and HBVM(18,2), 
respectively. Since the methods leave the [x, y)-plane invariant for the motion of the primaries and the 
z-axis invariant for the motion of the planetoid, we have just reported the motion of the primaries in the 
(x, y)-phase plane (left pictures) and the space-time diagram of the planetoid (right picture). 

We observe that, for the Gauss method, the orbits of the primaries are irregular in character so that 
the third body, after performing some oscillations around the origin, will eventually leave the system (see 
the right picture of Figure |3j). On the contrary (left picture of Figure [4]), the HBVM(18,2) generates a 
quite regular phase portrait. Due to the large stepsize h used, a sham rotation of the (x, y)-plane appears 
which, however, does not destroy the global symmetry of the dynamics, as testified by the bounded 
oscillations of the planetoid (right picture of Figure 0]) which look very similar to the reference ones in 
Figurc[TJ This aspect is also confirmed by the pictures in Figure[5j displaying the distance of the primaries 
as a function of the time. We see that the distance of the apocentres (corresponding to the maxima in 
the plots), as the two bodies wheel around the origin, are preserved by the HBVM(18,2) (right picture) 
while the same is not true for the Gauss method (left picture). 

3.2 The Henon-Heiles problem 

The Henon-Heiles equation originates from a problem in Celestial Mechanics describing the motion of 
a star under the action of a gravitational potential of a galaxy which is assumed time-independent and 
with an axis of symmetry (the z-axis) (see [llj and references therein). The main question related to 
this model was to state the existence of a third first integral, beside the total energy and the angular 
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Figure 6: Level curves of the potential U(qi,q2) of the Henon-Heiles problem (see (|14l) ). The origin O is 
a stable equilibrium point, whose domain of stability contains the equilateral triangle having as vertices 
the saddle points Pi, P%, and P3, provided that the total energy does not exceed the value \. Inside the 
triangle an orbit (qi(t), <Z2(i)) is traced whose total energy is close (but lower than) |. The trajectory 
gets very close to the sides of the triangle, which makes the problem of conserving the total energy in the 
numerical solution an important feature to avoid instability when a large stepsize is used. 

momentum^ By exploiting the symmetry of the system and the conservation of the angular momentum, 
Henon and Heiles reduced from three (cylindrical coordinates) to two (planar coordinates) the degrees 
of freedom, thus showing that the problem was equivalent to the study of the motion of a particle in a 
plane subject to an arbitrary potential U(q\, q%): 

H{q,p) = \{ P l+pl) + U{q 1: q 2 ). (13) 

Since U in (I13|) has no symmetry in general, we cannot consider the angular momentum as an invariant 
anymore, so that the only known first integral is the total energy represented by (fT3"]) itself, and the 
question is whether or not a second integral does exist. Henon and Heiles conducted a series of tests 
with the aim of giving a numerical evidence of the existence of such integral for moderate values of the 
energy H, and of the appearance of chaotic behavior when H{q,p) becomes larger than a critical value. 
In particular, for their experiments they choose 

ufata) = gfai + <&) + tin - 3«2> (14) 

which makes the Hamiltonian function a polynomial of degree three. 

When £/(<7i, 92) approaches the value i, the level curves of U tend to an equilateral triangle, whose 
vertices are saddle points of U (see FigureE]). This vertices have coordinates P\ = (0, 1), P2 = (— 2 > —5) 
andP 3 = (^,-i). 

We consider an initial point (q ,p ) such that q is inside the triangle U < ^ and H(q ,p ) < |: then 
the orbit originating from (q ,p Q ) will never abandon the triangle for any value of the time t. However, 
when H(q ,p ) is chosen very close to g, a numerical method which does not preserve exactly the total 

3 An analytical approach to the problem may be found in 1 10) . where the author finds out a formal expansion of the third 
invariant. 
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energy could cause the (numerical) orbit to jump outside the triangle and possibly to diverge to infinity. 
This aspect is further emphasized when a large stepsize of integration is used, as is usually required in 
the long time simulation of a dynamical system. 




Figure 7: The numerical trajectory in the (gi, (j2)-plane computed by the Gauss method of order four 
with stepsize h = 1. The stable character of the continuous orbit is not correctly reproduced by the 
numerical method: after a time t ~ 7000 the orbit escapes from the triangle (see the dots surrounded by 
small circles at the bottom right of the picture). 

We have integrated problem flT5|) in the time interval [0, 5 • 10 4 ] with stepsize h = 1 by using the Gauss 
method of order four (HBVM(2,2)) and the HBVM(4,2) method which assures an exact conservation of 
the total energy. 

Figures [7] and [8] show the numerical trajectories in the (qx, q^-plane as dots that eventually will 
densely fill the triangle. The orbit generated by the Gauss method is plotted up to time t ~ 7000, since 
it then escapes from the triangle, as highlighted by the three circles close to the saddle point P3. In 
fact, as Figure IH1 shows, the numerical Hamiltonian function associated with the Gauss method produces 
very irregular oscillations around the theoretical value (straight line) which eventually determine a loss 
of stability. 

On the contrary, all the 50000 dots of the numerical trajectory computed by the HBVM(4,2) method 
are visible in Figure [5] 

3.3 Computing the period annulus of a non- degenerate center of a polynomial Hamiltonian planar sys- 
tem. 

Non-degenerate centers^ of planar, in particular polynomial, Hamiltonian systems are extensively re- 
searched in the modern literature (see [91 [5] and references therein). The integration of such 
systems by means of HBVMs deserves a particular interest because, the degrees of freedom being one, 
the corresponding numerical solution is guaranteed to lie on the same level set H{q,p) = H(qq,pq) as the 
theoretical orbit. Furthermore, if this latter consists of a closed orbit surrounding an equilibrium point 
(center), the numerical solution will (in general) fill densely the corresponding closed level curve, thus 

4 We recall that a center is an equilibrium point which is surrounded by periodic orbits. It is non-degenerate if the 
linearized vector field at this point has non-zero eigenvalues. 
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Figure 8: The numerical trajectory in the (qi, <;2)-plane computed by the HBVM(4,2) method with 
stepsize h = 1. Since this method precisely conserves the total energy of the system, the orbit is entirely 
contained in the triangle at all times. 

reproducing the very same phase portrait associated with the original continuous problem. 

The region of marginal stability of a center Po, is called the period annulus of Pq and will be denoted 
by V: it is the largest punctured neighborhood of the center consisting of only periodic orbits. The 
function which associates to any periodic orbit in V its period is called the period function of the center. 
Such function has been being intensively studied for many years: its behavior relates to problems of 
isochronicityjf] monotonicity, bifurcation of its critical points, etc. 

The aim of the present example is to consider one such system and try to reproduce numerically, as 
best as possible, the set dV, that is the boundary of the period annulus V . Let H* < +00 be the value 
of the Hamiltonian function corresponding to any points on dP @ The Hamiltonian function we consider 
here is the fifth-degree polynomial 

H(p, q) = A(p) + B(p)q + C(p)q 2 + D(p)q\ (15) 

where 

A{P) = V 2 {\ + c 3 p + hp 2 + a 3 p 3 ), B(p) = p 2 (c 2 + b 2 p + a 2 p 2 ), 

C(p) = \ + cip + hp 2 + aip 3 , D{p) = c + b p + a p 2 , 

with (00,01,02,03) 7^ (0,0,0,0)0 Note that, since H(q,p) = \{p 2 + q 2 ) + h.o.t., we can assume P to be 
the origin O = (0, 0). 

The class of Hamiltonian systems defined by (|15|l has been proposed in 20 and [2"T]F1 Their main 
result was proving that the origin may not be an isochronous center [20 and, more specifically, that 
the period tends to infinity as H(qo,po) /*■ H*, (qo,Po) being the initial condition associated with the 
differential system. 

5 Namely, all the orbits surrounding the center Po share the same period. 

6 Here we assume that the center Po is non global: this is certainly true if H(q,p) is a polynomial of odd degree. 

7 Otherwise the degree of H(q,p) becomes lower than 5. 

8 The authors showed that, without loss of generality, the form (1 1 5 I t may be associated to any polynomial Hamiltonian 
system of degree four and admitting a non-degenerate center, via a suitable change of coordinates. 
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Figure 9: Hamiltonian function evaluated along the numerical solution of the Gauss and HBVM(4,2) 
methods. The irregular oscillations introduced by the Gauss method will cause the associated numerical 
solution to eventually leave the stability region centered at the origin. 

For our experiments, we have set the values of the coefficients {(Zj}, {6,}, and {ci} as follows: 



a = 0; 
bo = 0; 
c = 0; 



ai = 0; 
h = 1; 
Ci = 1; 



a 2 = 1 
b 2 = 
c 2 = 1 



a 3 = 0; 
h = 1; 
c 3 = 0. 



(16) 



In such a case, besides the origin Pq = (0,0), H(q,p) admits the following real equilibrium points (up to 
the machine precision): 

P l = (-6.879526475540134- 10"\ - 5.206527058470621 ■ 10" 1 ) — > saddle point; 
P 2 = (-1.179582379893681, 1.756351969248087) — > saddle point. 

Figure [TU] reports the shape of the level curves of ([T5l) (fTB"]) in a region enclosing Pq and P\ . We see that 
the limit closed orbit corresponding to dV is the one embracing Pq and having P\ as both w-limit point 
and a-limit poinlQ and, therefore, the value H* may be computed with precision as 

H* = H{Px) = 9.050199350868576 • 10~ 2 . (17) 

Now suppose we do not know the value H* in (fT7|) (it will be used as a reference value) and that 
we want to reproduce the orbit covering dV by simply picking initial points (qo,Po) further and further 
away from the origin, and checking whether the numerical solution remains bounded over a long timeF^I 
More precisely, we will locate the limit cycle by means of a dichotomic search, according to the following 
algorithm: 

step 1: find a point Q from which an orbit originates that does not embraces the critical point Pq (that is 

step 2: consider the segment joining Pq to Q: 

j(c) = (1-c)P + cPx, ce[0,l], 

and set cq = and c\ = 1; 

9 That is, limt—,±ca(q(t) , p(t)) = Pi for any choice of (go,Po) £ BP, 
10 Of course, we cannot assume (goiPo) = since Pi is an equilibrium point. 
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step 3: if c\ — Co < tol, STOP (tol is a specified tolerance); 

step 4: set c = CQ + Cl and solve numerically the Hamiltonian problem defined in (fl~5|) . considering 7(c) as 
initial condition, in the time interval [0, hN] where h > is the stepsize and N is a positive integer 
such that hN is large enough to give some information about the fate of the orbit originating from 
7(c). 

step 5: if the numerical solution eventually depart from Pq, set c\ = c, otherwise set cq — c, go to step 3; 

The point yo = (go, Po) = 7(c), where c is the value resulting after the execution of the above 
procedure, may be assumed as a point on dV within the specified tolerance tol. Detecting the limit cycle 
with high accuracy requires a huge number of simulations and therefore large run times, also taking into 
account the wide time intervals that must be used in order to inspect the asymptotic behavior of the 
numerical solutionis Consequently, it would be advisable to work with a relatively large stepsize h. We 
have set: 

h = 1, 0.5, N = 2500, 5000, tol = 2~ 52 (i.e., the value of eps in Matlab), Q = (0, 1), 
to cover the integration interval [0, 2500]. 
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Table 1: A point yo on the boundary of the period annulus V is computed by the Gauss and HBVM 
methods of orders 4, 6, 8 and 10 (corresponding to s = 2, 3, 4, 5 respectively). By their very nature, if 
used with a sufficient number of silent stages, HBVMs produce a numerical orbit that precisely lie on the 
same level set H(q,p) = H(qo,po) as the theoretical one, therefore we see that HBVMs can locate the 
point 2/0 with extreme precision, whatever the order and/or the stepsize used. On the contrary, Gauss 
methods produce a certain error that may be lowered by reducing the stepsize of integration h and/or 
by raising their order. 

Table Q] compares the results obtained by using the Gauss (HBVM(s,s)) and HBVM(fc,s) methods 
of orders 4, 6, 8 and 10 (therefore, since s — 2,3, 4, 5, we must choose, according to ©, k — 5,8, 10, 13, 
respectively, in order for the HBVM(fc,s) to exactly conserve the Hamiltonian function). We have denoted 
by y^'^ the point computed by the method HBVM(fc,s), and reported the error \H(yQ k ' s ^) — H*\/H* 
to estimate the accuracy with which each method computes the boundary of V . As was expected, the 
accuracy in detecting the right boundary of the period annulus by means of HBVMs is of the same order 
as the machine precision whatever the order and stepsize used (indeed, the value of y^'^ remains the 
same for all simulations). On the contrary, the Gauss methods produce a certain error which depends 

1 Actually, by virtue of their conservation properties, HBVMs do not need to be integrated over a long time, even though 
here we do that for comparison purposes. 
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both on the stepsize and on the order used: increasing the accuracy would require a suitable reduction of 
the stepsize and/or a grow-up of the order. Figure [TT] shows that even small oscillations of the numerical 
Hamiltonian function (left picture) could produce a noticeable irregularity of the numerical orbit in a 
neighborhood of the boundary of the period annulus (right picture). By their very nature, HBVMs 
succeed in detecting the set dV with an accuracy of the same order as the machine precision: the error 
in the Hamiltonian function is negligible (left picture) and the numerical orbit correctly passes through 
the saddle point Pi. 
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Figure 10: Level curves of the Hamiltonian (jT5j) in a region that embraces the center point Po and the 
saddle point Pi . Each level curve, corresponding to an orbit of the associated Hamiltonian system, is 
labeled by a number that indicates its elevation. 
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